Корреляция

Необходимые библиотеки:

library(dplyr)
library(ggplot2)
library(PerformanceAnalytics)
library(RColorBrewer)
library(corrplot)
library(Hmisc)
library(heatmaply)

Данные

Данные содержат некоторые индивидуальные характеристики студентов и результаты их успеваемости по курсу.

Переменные

• user_id (int) - идентификатор пользователя

• gender (factor) - пол

  male	female	NA
  мужской	женский	не указан

• age (int) - возраст
• grade (int) - оценка за курс (диапазон принимаемых значений: от 0 до 100)
• number_lectures (int) - количество лекций, просмотренных студентом
• number_supplement_materials (int) - количество дополнительных материалов, просмотренных студентом
• number_forum_posts (int) - количество сообщений на форуме, написанных студентом

• number_pauses (int) - количество пауз, сделанных студентом в видеолекциях

---

Считываем данные

data <- read.csv("data.csv")
head(data)

Проверяем все ли типы переменных соответствуют нужным

data$user_id <- as.integer(data$user_id)
data$gender <- as.factor(data$gender)
data$age <- as.integer(data$age)
data$grade <- as.integer(data$grade)
data$number_lectures <- as.integer(data$number_lectures)
data$number_supplement_materials <- as.integer(data$number_supplement_materials)
data$number_forum_posts <- as.integer(data$number_forum_posts)
data$number_pauses <- as.integer(data$number_pauses)

Посмотрим на данные

Пропущенные значения

Посмотрим есть ли в данных пропущенные значения. Если есть, удалим такие наблюдения (строки).

head(is.na(data))

##      user_id gender   age grade number_lectures number_supplement_materials
## [1,]   FALSE  FALSE FALSE FALSE           FALSE                       FALSE
## [2,]   FALSE  FALSE FALSE FALSE           FALSE                       FALSE
## [3,]   FALSE  FALSE FALSE FALSE           FALSE                       FALSE
## [4,]   FALSE  FALSE FALSE FALSE           FALSE                       FALSE
## [5,]   FALSE  FALSE FALSE FALSE           FALSE                       FALSE
## [6,]   FALSE  FALSE FALSE FALSE           FALSE                       FALSE
##      number_forum_posts number_pauses
## [1,]              FALSE         FALSE
## [2,]              FALSE         FALSE
## [3,]              FALSE         FALSE
## [4,]              FALSE         FALSE
## [5,]              FALSE         FALSE
## [6,]              FALSE         FALSE

sum(is.na(data))

## [1] 37

data_with_NA <- data
data <- na.omit(data)

Выбросы

Посмотрим есть ли выбросы в данных.

boxplot(data)

boxplot(data$grade)

Удалим выбросы в переменной grade

outliers <- boxplot(data$grade, plot=FALSE)$out
data_with_outliers <- data
data <- data[-which(data$grade %in% outliers),]
boxplot(data$grade)

Посмотрим выбросы в переменной number_lectures

boxplot(data$number_lectures)

Не будем удалять выбросы, так как ошибочных или слишком сильно отклоняющихся значений в переменной number_lectures нет

Допущения

Необходимо помнить два важных допущения для построения корреляции Пирсона:

• Обе переменные распределены нормально
• Связь между переменными линейна

Посмотрим, как распределены некоторые переменные с помощью диаграммы частот. И проверим нормальность распределения с помощью теста Шапиро-Уилка.

hist(data$grade, breaks = 20, col = "lightblue")

shapiro.test(data$grade)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$grade
## W = 0.99626, p-value = 0.2927

Диаграмма показывает, что условие нормальности распределения выполняется. p-value = 0.2927, т.е. p-value > 0.05. Соответственно, гипотеза о нормальности распределения не отвергается.

hist(data$number_lectures, breaks = 20, col = "lightblue")

shapiro.test(data$number_lectures)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$number_lectures
## W = 0.99477, p-value = 0.08828

Диаграмма показывает, что условие нормальности распределения выполняется. p-value = 0.08828, т.е. p-value < 0.05. Соответственно, гипотеза о нормальности распределения не отвергается.

hist(data$age, breaks = 20, col = "lightblue")

shapiro.test(data$age)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$age
## W = 0.94899, p-value = 4.302e-12

Диаграмма показывает, что условие нормальности распределения не выполняется. p-value = 4.302e-12, т.е. p-value < 0.05. Соответственно, гипотеза о нормальности распределения отвергается.

Построим корреляцию на двух переменных

Визуализация

Диаграмма рассеяния (scatterplot)

plot(data$number_lectures, 
     data$grade, 
     type = "p",
     main = "Соотношение числа просмотренных лекций с оценкой за курс",
     xlab = "Количество просмотренных лекций",
     ylab = "Оценка за курс",
     col = "steelblue",
     pch = 19)

ИЛИ

ggplot(data, aes(x = number_lectures, y = grade)) + 
    geom_point(color="#69b3a2") +
  xlab("Количество просмотренных лекций") +
  ylab("Оценка за курс")

Добавляем прямую линию (или линию регрессии. О ней более подробно будет с следующей лекции)

ggplot(data, aes(x = number_lectures, y = grade)) +
  geom_point(color="#69b3a2") +
  xlab("Количество просмотренных лекций") +
  ylab("Оценка за курс") +
  geom_smooth(method = lm, color = "red", se = FALSE)

## `geom_smooth()` using formula 'y ~ x'

Расчет коэффициента корреляции

В R коэффициент корреляции можно легко рассчитать с помощью функциий cor() и cor.test(). Различие между данными функциями состоит в том, что cor() позволяет расчитать только коэффициент корреляции, а cor.test() помимо расчета самого коэффициента также выводит статистическую значимость коэффициента.

1. Рассчитам коэффициент корреляции между переменными grade и number_lectures.

Мы используем корреляцию Пирсона, так как оба допущения выполняются: обе переменные распределены нормально и связь между переменными линейна

cor(data$grade, data$number_lectures)

## [1] 0.6757337

Коэффициент корреляции Пирсона равен 0.68, что указывает на положительную корреляцию между количеством проссмотренных студентом лекций и его оценкой за курс.

cor.test(data$grade, data$number_lectures)

## 
##  Pearson's product-moment correlation
## 
## data:  data$grade and data$number_lectures
## t = 20.436, df = 497, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6250288 0.7207604
## sample estimates:
##       cor 
## 0.6757337

Во втором случае коэффициент корреляции Пирсона тоже равен 0.68, что указывает на положительную корреляцию. Также p-value сотавляет 2.2e-16, что меньше 0.05. Соответственно, гипотеза об отсутствии корреляции отвергается. Значит корреляция между количеством проссмотренных студентом лекций и его оценкой за курс есть.

2. Рассчитам коэффициент корреляции между переменными grade и age.

Поскольку в нашем случае не выполняется условие нормальности распределения для обеих переменных, необходимо использовать непараметрический коэффициент корреляции Спирмена. Для вычисления коэффициента корреляции Спирмена в R также используются функции cor() и cor.test(). Однако, по умолчанию данные функции расчитывают коэффициент корреляции Пирсона. Чтобы расчитать коэффициент корреляции Спирмена необходимо установить значение “spearman” в аргументе method:

cor(data$grade, data$age, method = "spearman")

## [1] 0.03545132

Коэффициент корреляции Спирмена равен 0.04, что указывает на очень слабую (почти нулевую) корреляцию между возрастом студента и его оценкой за курс. В таком случае можно скачать, что корреляция отсутствует.

cor.test(data$grade, data$age, method = "spearman")

## Warning in cor.test.default(data$grade, data$age, method = "spearman"): Cannot
## compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$grade and data$age
## S = 19974356, p-value = 0.4294
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.03545132

Во втором случае коэффициент корреляции Спирмена тоже равен 0.04. Также p-value сотавляет 0.4294, что больше 0.05. Соответственно, гипотеза об отсутствии корреляции не отвергается. Значит корреляции между возрастом студентов и их оценкой за курс нет.

Корреляционная матрица

Визуализация

1. Базовая матрица диаграмм рассеяния (Scatterplot Matrix)

pairs( ~ age + grade + number_lectures + number_pauses, data = data,
       main="Scatterplot Matrix")

Построение матрицы

Важно помнить про допущения корреляции. Проверим распределение переменных age, number_supplement_materials, number_forum_posts и number_pauses с помощью теста Шапиро-Уилка.

shapiro.test(data$age)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$age
## W = 0.94899, p-value = 4.302e-12

shapiro.test(data$number_supplement_materials)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$number_supplement_materials
## W = 0.99347, p-value = 0.02953

shapiro.test(data$number_forum_posts)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$number_forum_posts
## W = 0.98541, p-value = 6.705e-05

shapiro.test(data$number_pauses)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$number_pauses
## W = 0.96991, p-value = 1.352e-08

Как мы видим, во всех случаях p-value < 0.05. Соответственно, гипотеза о нормальности распределения отвергается. Значит при построении корреляционной матрицы мы далжны указать, что method = “spearman”.

data_1 <- data[,c(3,6,7,8)]
cor(data_1, method = "spearman")

##                                    age number_supplement_materials
## age                         1.00000000                  0.07663774
## number_supplement_materials 0.07663774                  1.00000000
## number_forum_posts          0.01857328                  0.51466494
## number_pauses               0.05511923                  0.17085916
##                             number_forum_posts number_pauses
## age                                 0.01857328    0.05511923
## number_supplement_materials         0.51466494    0.17085916
## number_forum_posts                  1.00000000    0.21436328
## number_pauses                       0.21436328    1.00000000

Другие варианты визуализации

1. heatmap (или тепловая карта) с помощью функции heatmaply_cor()

heatmaply_cor(
  cor(data_1, method = "spearman"),
  Colv=NA, Rowv=NA, symm = TRUE
)

## Warning in doTryCatch(return(expr), name, parentenv, handler): не могу загрузить разделяемый объект '/Library/Frameworks/R.framework/Resources/modules//R_X11.so':
##   dlopen(/Library/Frameworks/R.framework/Resources/modules//R_X11.so, 6): Library not loaded: /opt/X11/lib/libSM.6.dylib
##   Referenced from: /Library/Frameworks/R.framework/Versions/3.6/Resources/modules/R_X11.so
##   Reason: image not found

• Аргумент Colv = NA означает, что мы хотим получить простой heatmap (без веток - без опрделения порядка переменных и кластеризации)
• Аргумент Rowv = NA - мы хотим получить простой heatmap (без веток - без опрделения порядка переменных и кластеризации)
• Аргумент symm = TRUE - логическое указание, следует ли рассматривать ось X симметрично

2. corrplot

cor_corplot <- cor(data_1, method = "spearman")
corrplot(cor_corplot, method = "circle")

corrplot(cor_corplot, method = "number")

3. chart.Correlation

chart.Correlation(data_1, histogram=TRUE, pch=19)

Домашняя работа

Попробуйте построить рассчитать корреляцию между переменными

1. Выберете 2 заинтересовавшие Вас переменные

2. Проверьте для них допущения

3. Постройте диаграмму рассеяния и рассчитайте коэффициент корреляции

4. Проинтерпретируйте полученные результаты

5. Помимо этих двух переменных, отберите еще 1-2 заинтересовавших Вас переменных в новый датафрейм

6. Постройте корреляционную любым способом

7. Визуализируйте корреляционную матрицу любыми 2-мя способами